/* v5d.c */

/* Vis5D version 5.2 */

/*
Vis5D system for visualizing five dimensional gridded data sets
Copyright (C) 1990 - 2000 Bill Hibbard, Johan Kellum, Brian Paul,
Dave Santek, and Andre Battaiola.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 1, or (at your option)
any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/


/* this should be updated when the file version changes */
#define FILE_VERSION "4.3"



/*
 * New grid file format for VIS-5D:
 *
 * The header is a list of tagged items.  Each item has 3 parts:
 *    1. A tag which is a 4-byte integer identifying the type of item.
 *    2. A 4-byte integer indicating how many bytes of data follow.
 *    3. The binary data.
 *
 * If we need to add new information to a file header we just create a
 * new tag and add the code to read/write the information.
 *
 * If we're reading a header and find an unknown tag, we can use the
 * length field to skip ahead to the next tag.  Therefore, the file
 * format is forward (and backward) compatible.
 *
 * Grid data is stored as either:
 *     1-byte unsigned integers  (255=missing)
 *     2-byte unsigned integers  (65535=missing)
 *     4-byte IEEE floats        ( >1.0e30 = missing) 
 *
 * All numeric values are stored in big endian order.  All floating point
 * values are in IEEE format.
 */



/*
 * Updates:
 *
 * April 13, 1995, brianp
 *   finished Cray support for 2-byte and 4-byte compress modes
 */




#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "binio.h"
#include "v5d.h"
#include "vis5d.h"
#ifndef SEEK_SET
#  define SEEK_SET 0
#endif
#ifndef SEEK_CUR
#  define SEEK_CUR 1
#endif
#ifndef SEEK_END
#  define SEEK_END 2
#endif



/*
 * Currently defined tags:
 * Note:  the notation a[i] doesn't mean a is an array of i elements,
 * rather it just refers to the ith element of a[].
 *
 * Tags marked as PHASED OUT should be readable but are no longer written.
 * Old tag numbers can't be reused!
 * 
 */


/*      TAG NAME        VALUE       DATA (comments)                     */
/*----------------------------------------------------------------------*/
#define TAG_ID          0x5635440a  /* hex encoding of "V5D\n"          */

/* general stuff 1000+ */
#define TAG_VERSION     1000        /* char*10 FileVersion              */
#define TAG_NUMTIMES    1001        /* int*4 NumTimes                   */
#define TAG_NUMVARS     1002        /* int*4 NumVars                    */
#define TAG_VARNAME     1003        /* int*4 var; char*10 VarName[var]  */

#define TAG_NR          1004        /* int*4 Nr                         */ 
#define TAG_NC          1005        /* int*4 Nc                         */
#define TAG_NL          1006        /* int*4 Nl  (Nl for all vars)      */
#define TAG_NL_VAR      1007        /* int*4 var; int*4 Nl[var]         */
#define TAG_LOWLEV_VAR  1008        /* int*4 var; int*4 LowLev[var]     */

#define TAG_TIME        1010        /* int*4 t;  int*4 TimeStamp[t]     */
#define TAG_DATE        1011        /* int*4 t;  int*4 DateStamp[t]     */

#define TAG_MINVAL      1012        /* int*4 var;  real*4 MinVal[var]   */
#define TAG_MAXVAL      1013        /* int*4 var;  real*4 MaxVal[var]   */

#define TAG_COMPRESS    1014        /* int*4 CompressMode; (#bytes/grid)*/

#define TAG_UNITS       1015        /* int *4 var; char*20 Units[var]   */

/* vertical coordinate system 2000+ */
#define TAG_VERTICAL_SYSTEM 2000    /* int*4 VerticalSystem             */
#define TAG_VERT_ARGS    2100       /* int*4 n;  real*4 VertArgs[0..n-1]*/

#define TAG_BOTTOMBOUND  2001       /* real*4 BottomBound     (PHASED OUT)  */
#define TAG_LEVINC       2002       /* real*4 LevInc      (PHASED OUT)      */
#define TAG_HEIGHT       2003    /* int*4 l;  real*4 Height[l] (PHASED OUT) */


/* projection 3000+ */
#define TAG_PROJECTION   3000       /* int*4 projection:                    */
                                    /*   0 = generic linear                 */
                                    /*   1 = cylindrical equidistant        */
                                    /*   2 = Lambert conformal/Polar Stereo */
                                    /*   3 = rotated equidistant            */
#define TAG_PROJ_ARGS    3100       /* int *4 n;  real*4 ProjArgs[0..n-1]   */

#define TAG_NORTHBOUND   3001       /* real*4 NorthBound       (PHASED OUT) */
#define TAG_WESTBOUND    3002       /* real*4 WestBound        (PHASED OUT) */
#define TAG_ROWINC       3003       /* real*4 RowInc           (PHASED OUT) */
#define TAG_COLINC       3004       /* real*4 ColInc           (PHASED OUT) */
#define TAG_LAT1         3005       /* real*4 Lat1             (PHASED OUT) */
#define TAG_LAT2         3006       /* real*4 Lat2             (PHASED OUT) */
#define TAG_POLE_ROW     3007       /* real*4 PoleRow          (PHASED OUT) */
#define TAG_POLE_COL     3008       /* real*4 PoleCol          (PHASED OUT) */
#define TAG_CENTLON      3009       /* real*4 CentralLon       (PHASED OUT) */
#define TAG_CENTLAT      3010       /* real*4 CentralLat       (PHASED OUT) */
#define TAG_CENTROW      3011       /* real*4 CentralRow       (PHASED OUT) */
#define TAG_CENTCOL      3012       /* real*4 CentralCol       (PHASED OUT) */
#define TAG_ROTATION     3013       /* real*4 Rotation         (PHASED OUT) */
#define TAG_ROWINCKM     3014       /* real*4 RowIncKm         (PHASED OUT) */
#define TAG_COLINCKM     3015       /* real*4 ColIncKm         (PHASED OUT) */


#define TAG_END                9999






/**********************************************************************/
/*****                  Miscellaneous Functions                   *****/
/**********************************************************************/


float pressure_to_height(float pressure)
{
  return (float) DEFAULT_LOG_EXP * log((double) pressure / DEFAULT_LOG_SCALE);
}

float height_to_pressure(float height)
{
  return (float) DEFAULT_LOG_SCALE * exp((double) height / DEFAULT_LOG_EXP);
}


/*
 * Return current file position.
 * Input:  f - file descriptor
 */
static off_t ltell( int f )
{
   return lseek( f, 0, SEEK_CUR );
}


/*
 * Copy up to maxlen characters from src to dst stopping upon whitespace
 * in src.  Terminate dst with null character.
 * Return:  length of dst.
 */
static int copy_string2( char *dst, const char *src, int maxlen )
{
   int i;

   for (i=0;i<maxlen;i++) dst[i] = src[i];
   for (i=maxlen-1; i>=0; i--) {
     if (dst[i]==' ' || i==maxlen-1) dst[i] = 0;
     else break;
   }
   return strlen(dst);
}



/*
 * Copy up to maxlen characters from src to dst stopping upon whitespace
 * in src.  Terminate dst with null character.
 * Return:  length of dst.
 */
static int copy_string( char *dst, const char *src, int maxlen )
{
   int i;

   for (i=0;i<maxlen;i++) {
      if (src[i]==' ' || i==maxlen-1) {
         dst[i] = 0;
         break;
      }
      else {
         dst[i] = src[i];
      }
   }
   return i;
}



/*
 * Convert a date from YYDDD format to days since Jan 1, 1900.
 */
int v5dYYDDDtoDays( int yyddd )
{
   int iy, id, idays;

   iy = yyddd / 1000;
   id = yyddd - 1000*iy;
   if (iy < 50) iy += 100; /* WLH 31 July 96 << 31 Dec 99 */
   idays = 365*iy + (iy-1)/4 + id;

   return idays;
}


/*
 * Convert a time from HHMMSS format to seconds since midnight.
 */
int v5dHHMMSStoSeconds( int hhmmss )
{
   int h, m, s;

   h = hhmmss / 10000;
   m = (hhmmss / 100) % 100;
   s = hhmmss % 100;

   return s + m*60 + h*60*60;
}



/*
 * Convert a day since Jan 1, 1900 to YYDDD format.
 */
int v5dDaysToYYDDD( int days )
{
   int iy, id, iyyddd;

   iy = (4*days)/1461;
   id = days-(365*iy+(iy-1)/4);
   if (iy > 99) iy = iy - 100; /* WLH 31 July 96 << 31 Dec 99 */
   /* iy = iy + 1900; is the right way to fix this, but requires
      changing all places where dates are printed - procrastinate */
   iyyddd = iy*1000+id;

   return iyyddd;
}


/*
 * Convert a time in seconds since midnight to HHMMSS format.
 */
int v5dSecondsToHHMMSS( int seconds )
{
   int hh, mm, ss;

   hh = seconds / (60*60);
   mm = (seconds / 60) % 60;
   ss = seconds % 60;
   return hh*10000 + mm * 100 + ss;
}




void v5dPrintStruct( const v5dstruct *v )
{
   static char day[7][10] = { "Sunday", "Monday", "Tuesday", "Wednesday",
                              "Thursday", "Friday", "Saturday" };
   int time, var, i;
   int maxnl;

   maxnl = 0;
   for (var=0;var<v->NumVars;var++) {
      if (v->Nl[var]+v->LowLev[var]>maxnl) {
         maxnl = v->Nl[var]+v->LowLev[var];
      }
   }

   if (v->FileFormat==0) {
      if (v->FileVersion[0] == 0) {
        printf("File format: v5d  version: (4.0 or 4.1)\n");
      }
      else {
        printf("File format: v5d  version: %s\n", v->FileVersion);
      }
   }
   else {
      printf("File format: comp5d  (VIS-5D 3.3 or older)\n");
   }

   if (v->CompressMode==1) {
      printf("Compression:  1 byte per gridpoint.\n");
   }
   else {
      printf("Compression:  %d bytes per gridpoint.\n", v->CompressMode);
   }
   printf("header size=%d\n", v->FirstGridPos);
   printf("sizeof(v5dstruct)=%d\n", sizeof(v5dstruct) );
   printf("\n");

   printf("NumVars = %d\n", v->NumVars );

   printf("Var  Name       Units      Rows  Cols  Levels LowLev  MinVal       MaxVal\n");
   for (var=0;var<v->NumVars;var++) {
      printf("%3d  %-10s %-10s %3d   %3d   %3d    %3d",
             var+1, v->VarName[var], v->Units[var],
             v->Nr, v->Nc, v->Nl[var], v->LowLev[var] );
      if (v->MinVal[var] > v->MaxVal[var]) {
         printf("     MISSING      MISSING\n");
      }
      else {
         printf("     %-12g %-12g\n", v->MinVal[var], v->MaxVal[var] );
      }
   }

   printf("\n");

   printf("NumTimes = %d\n", v->NumTimes );
   printf("Step    Date(YYDDD)    Time(HH:MM:SS)   Day\n");
   for (time=0;time<v->NumTimes;time++) {
      int i = v->TimeStamp[time];
      printf("%3d        %05d       %5d:%02d:%02d     %s\n",
             time+1,
             v->DateStamp[time],
             i/10000, (i/100)%100, i%100,
             day[ v5dYYDDDtoDays(v->DateStamp[time]) % 7 ]);
   }
   printf("\n");

   switch (v->VerticalSystem) {
      case 0:
         printf("Generic linear vertical coordinate system:\n");
         printf("\tBottom Bound: %f\n", v->VertArgs[0] );
         printf("\tIncrement between levels:  %f\n", v->VertArgs[1] );
         break;
      case 1:
         printf("Equally spaced levels in km:\n");
         printf("\tBottom Bound: %f\n", v->VertArgs[0] );
         printf("\tIncrement: %f\n", v->VertArgs[1] );
         break;
      case 2:
         printf("Unequally spaced levels in km:\n");
         printf("Level\tHeight(km)\n");
         for (i=0;i<maxnl;i++) {
            printf("%3d     %10.3f\n", i+1, v->VertArgs[i] );
         }
         break;
      case 3:
         printf("Unequally spaced levels in mb:\n");
         printf("Level\tPressure(mb)\n");
         for (i=0;i<maxnl;i++) {
            printf("%3d     %10.3f\n", i+1, height_to_pressure(v->VertArgs[i]) );
         }
         break;
      default:
         printf("Bad VerticalSystem value: %d\n", v->VerticalSystem );
   }
   printf("\n");

   switch (v->Projection) {
      case 0:
         printf("Generic linear projection:\n");
         printf("\tNorth Boundary: %f\n", v->ProjArgs[0] );
         printf("\tWest Boundary: %f\n", v->ProjArgs[1] );
         printf("\tRow Increment: %f\n", v->ProjArgs[2] );
         printf("\tColumn Increment: %f\n", v->ProjArgs[3] );
         break;
      case 1:
         printf("Cylindrical Equidistant projection:\n");
         printf("\tNorth Boundary: %f degrees\n", v->ProjArgs[0] );
         printf("\tWest Boundary: %f degrees\n", v->ProjArgs[1] );
         printf("\tRow Increment: %f degrees\n", v->ProjArgs[2] );
         printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
/*
         printf("\tSouth Boundary: %f degrees\n",
                v->NorthBound - v->RowInc * (v->Nr-1) );
         printf("\tEast Boundary: %f degrees\n",
                v->WestBound - v->ColInc * (v->Nc-1) );
*/
         break;
      case 2:
         printf("Lambert Conformal projection:\n");
         printf("\tStandard Latitude 1: %f\n", v->ProjArgs[0] );
         printf("\tStandard Latitude 2: %f\n", v->ProjArgs[1] );
         printf("\tNorth/South Pole Row: %f\n", v->ProjArgs[2] );
         printf("\tNorth/South Pole Column: %f\n", v->ProjArgs[3] );
         printf("\tCentral Longitude: %f\n", v->ProjArgs[4] );
         printf("\tColumn Increment: %f km\n", v->ProjArgs[5] );
         break;
      case 3:
         printf("Stereographic:\n");
         printf("\tCenter Latitude: %f\n", v->ProjArgs[0] );
         printf("\tCenter Longitude: %f\n", v->ProjArgs[1] );
         printf("\tCenter Row: %f\n", v->ProjArgs[2] );
         printf("\tCenter Column: %f\n", v->ProjArgs[3] );
         printf("\tColumn Spacing: %f\n", v->ProjArgs[4] );
         break;
      case 4:
         /* WLH 4-21-95 */
         printf("Rotated equidistant projection:\n");
         printf("\tLatitude of grid(0,0): %f\n", v->ProjArgs[0] );
         printf("\tLongitude of grid(0,0): %f\n", v->ProjArgs[1] );
         printf("\tRow Increment: %f degress\n", v->ProjArgs[2] );
         printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
         printf("\tCenter Latitude: %f\n", v->ProjArgs[4] );
         printf("\tCenter Longitude: %f\n", v->ProjArgs[5] );
         printf("\tRotation: %f degrees\n", v->ProjArgs[6] );
         break;
      case 5:
         printf("Mercator:\n");
         printf("\tCenter Latitude: %f\n", v->ProjArgs[0] );
         printf("\tCenter Longitude: %f\n", v->ProjArgs[1] );
         printf("\tRow Increment: %f Kilometers\n", v->ProjArgs[2] );
         printf("\tCol Increment: %f Kilometers\n", v->ProjArgs[3] );
         break;
      default:
         printf("Bad projection number: %d\n", v->Projection );
   }
}



/*
 * Compute the location of a compressed grid within a file.
 * Input:  v - pointer to v5dstruct describing the file header.
 *         time, var - which timestep and variable.
 * Return:  file offset in bytes
 */
static int grid_position( const v5dstruct *v, int time, int var )
{
   int pos, i;

   assert( time >= 0 );
   assert( var >= 0 );
   assert( time < v->NumTimes );
   assert( var < v->NumVars );

   pos = v->FirstGridPos + time * v->SumGridSizes;
   for (i=0;i<var;i++) {
      pos += v->GridSize[i];
   }

   return pos;
}



/*
 * Compute the ga and gb (de)compression values for a grid.
 * Input:  nr, nc, nl - size of grid
 *         data - the grid data
 *         ga, gb - arrays to store results.
 *         minval, maxval - pointer to floats to return min, max values
 *         compressmode - 1, 2 or 4 bytes per grid point
 * Output:  ga, gb - the (de)compression values
 *          minval, maxval - the min and max grid values
 * Side effect:  the MinVal[var] and MaxVal[var] fields in g may be
 *               updated with new values.
 */
static void compute_ga_gb( int nr, int nc, int nl, 
                           const float data[], int compressmode,
                           float ga[], float gb[],
                           float *minval, float *maxval )
{
#ifdef SIMPLE_COMPRESSION
   /*
    * Compute ga, gb values for whole grid.
    */
   int i, lev, allmissing, num;
   float min, max, a, b;

   min = 1.0e30;
   max = -1.0e30;
   num = nr * nc * nl;
   allmissing = 1;
   for (i=0;i<num;i++) {
      if (!IS_MISSING(data[i])) {
         if (data[i]<min)  min = data[i];
         if (data[i]>max)  max = data[i];
         allmissing = 0;
      }
   }
   if (allmissing) {
      a = 1.0;
      b = 0.0;
   }
   else {
      a = (max-min) / 254.0;
      b = min;
   }

   /* return results */
   for (i=0;i<nl;i++) {
      ga[i] = a;
      gb[i] = b;
   }

   *minval = min;
   *maxval = max;
#else
   /*
    * Compress grid on level-by-level basis.
    */
#  define SMALLVALUE -1.0e30
#  define BIGVALUE 1.0e30
#  define ABS(x)   ( ((x) < 0.0) ? -(x) : (x) )
   float gridmin, gridmax;
   float levmin[MAXLEVELS], levmax[MAXLEVELS];
   float d[MAXLEVELS], dmax;
   float ival, mval;
   int j, k, lev, nrnc;

   nrnc = nr * nc;

   /* find min and max for each layer and the whole grid */
   gridmin = BIGVALUE;
   gridmax = SMALLVALUE;
   j = 0;


   for (lev=0;lev<nl;lev++) {
      float ave, var;
      float min, max;
      min = BIGVALUE;
      max = SMALLVALUE;
      ave = 0.0;
      var = 0.0;
      for (k=0;k<nrnc;k++) {
         if (!IS_MISSING(data[j]) && data[j]<min)
            min = data[j];
         if (!IS_MISSING(data[j]) && data[j]>max)
            max = data[j];
         j++;
      }

      if (min<gridmin)
        gridmin = min;
      if (max>gridmax)
        gridmax = max;
      levmin[lev] = min;
      levmax[lev] = max;
   }

/* WLH 2-2-95 */
#ifdef KLUDGE
   /* if the grid minimum is within delt of 0.0, fudge all values */
   /* within delt of 0.0 to delt, and recalculate mins and maxes */
   {
      float delt;
      int nrncnl = nrnc * nl;

      delt = (gridmax - gridmin)/100000.0;
      if ( ABS(gridmin) < delt && gridmin!=0.0 && compressmode != 4 ) {
         float min, max;
         for (j=0; j<nrncnl; j++) {
            if (!IS_MISSING(data[j]) && data[j]<delt)
              data[j] = delt;
         }
         /* re-calculate min and max for each layer and the whole grid */
         gridmin = delt;
         for (lev=0;lev<nl;lev++) {
            if (ABS(levmin[lev]) < delt)
              levmin[lev] = delt;
            if (ABS(levmax[lev]) < delt)
              levmax[lev] = delt;
         }
      }
   }
#endif

   /* find d[lev] and dmax = MAX( d[0], d[1], ... d[nl-1] ) */
   dmax = 0.0;
   for (lev=0;lev<nl;lev++) {
      if (levmin[lev]>=BIGVALUE && levmax[lev]<=SMALLVALUE) {
         /* all values in the layer are MISSING */
         d[lev] = 0.0;
      }
      else {
         d[lev] = levmax[lev]-levmin[lev];
      }
      if (d[lev]>dmax)
         dmax = d[lev];
   }

   /*** Compute ga (scale) and gb (bias) for each grid level */
   if (dmax==0.0) {
      /*** Special cases ***/
      if (gridmin==gridmax) {
         /*** whole grid is of same value ***/
         for (lev=0; lev<nl; lev++) {
            ga[lev] = gridmin;
            gb[lev] = 0.0;
         }
      }
      else {
         /*** every layer is of a single value ***/
         for (lev=0; lev<nl; lev++) {
            ga[lev] = levmin[lev];
            gb[lev] = 0.0;
         }
      }
   }
   else {
      /*** Normal cases ***/
      if (compressmode == 1) {
#define ORIGINAL
#ifdef ORIGINAL
         ival = dmax / 254.0;
         mval = gridmin;

         for (lev=0; lev<nl; lev++) {
            ga[lev] = ival;
            gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival ); 
         }
#else
         for (lev=0; lev<nl; lev++) {
            if (d[lev]==0.0) {
               ival = 1.0;
            }
            else {
               ival = d[lev] / 254.0;
            }
            ga[lev] = ival;
            gb[lev] = levmin[lev];
         }
#endif
      }
      else if (compressmode == 2) {
         ival = dmax / 65534.0;
         mval = gridmin;

         for (lev=0; lev<nl; lev++) {
            ga[lev] = ival;
            gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival );
         }
      }
      else {
         assert( compressmode==4 );
         for (lev=0; lev<nl; lev++) {
            ga[lev] = 1.0;
            gb[lev] = 0.0;
         }
      }
   }

   /* update min, max values */
   *minval = gridmin;
   *maxval = gridmax;
#endif
}




/*
 * Compress a 3-D grid from floats to 1-byte unsigned integers.
 * Input: nr, nc, nl - size of grid
 *        compressmode - 1, 2 or 4 bytes per grid point
 *        data - array of [nr*nc*nl] floats
 *        compdata - pointer to array of [nr*nc*nl*compressmode] bytes
 *                   to put results into.
 *        ga, gb - pointer to arrays to put ga and gb decompression values
 *        minval, maxval - pointers to float to return min & max values
 * Output:  compdata - the compressed grid data
 *          ga, gb - the decompression values
 *          minval, maxval - the min and max grid values
 */
void v5dCompressGrid( int nr, int nc, int nl, int compressmode,
                      const float data[],
                      void *compdata, float ga[], float gb[],
                      float *minval, float *maxval )
{
   int nrnc = nr * nc;
   int nrncnl = nr * nc * nl;
   V5Dubyte *compdata1 = (V5Dubyte *) compdata;
   V5Dushort *compdata2 = (V5Dushort *) compdata;

   /* compute ga, gb values */
   compute_ga_gb( nr, nc, nl, data, compressmode, ga, gb, minval, maxval );

   /* compress the data */
   if (compressmode==1) {
      int i, lev, p;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float one_over_a, b;
/* WLH 5 Nov 98
         b = gb[lev] - 0.0001;
*/
         /* WLH 5 Nov 98 */
         b = gb[lev];
                                /* subtract an epsilon so the int((d-b)/a) */
                                /* expr below doesn't get mis-truncated. */
         if (ga[lev]==0.0) {
            one_over_a = 1.0;
         }
         else {
            one_over_a = 1.0 / ga[lev];
         }
         for (i=0;i<nrnc;i++,p++) {
            if (IS_MISSING(data[p])) {
               compdata1[p] = 255;
            }
            else {
/* MJK 1.19.99
               compdata1[p] = (V5Dubyte) (int) ((data[p]-b) * one_over_a);
*/
               compdata1[p] = (V5Dubyte) rint((data[p]-b) * one_over_a);
               if (compdata1[p] >= 255){
                  compdata1[p] = (V5Dubyte) (int) (255.0 - .0001);
               }
            }
         }
      }
   }

   else if (compressmode == 2) {
      int i, lev, p;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float one_over_a, b;
/* WLH 5 Nov 98
         b = gb[lev] - 0.0001;
*/
         /* WLH 5 Nov 98 */
         b = gb[lev];

         if (ga[lev]==0.0) {
            one_over_a = 1.0;
         }
         else {
            one_over_a = 1.0 / ga[lev];
         }
#ifdef _CRAY
         /* this is tricky because sizeof(V5Dushort)==8, not 2 */
         for (i=0;i<nrnc;i++,p++) {
            V5Dushort compvalue;
            if (IS_MISSING(data[p])) {
               compvalue = 65535;
            }
            else {
/* MJK 3.2.99
               compvalue = (V5Dushort) (int) ((data[p]-b) * one_over_a);
*/
               compvalue = (V5Dushort) rint((data[p]-b) * one_over_a);
            }
            compdata1[p*2+0] = compvalue >> 8;     /* upper byte */
            compdata1[p*2+1] = compvalue & 0xffu;  /* lower byte */
         }
#else
         for (i=0;i<nrnc;i++,p++) {
            if (IS_MISSING(data[p])) {
               compdata2[p] = 65535;
            }
            else {
               compdata2[p] = (V5Dushort) rint((data[p]-b) * one_over_a);

/*
               compdata2[p] = (V5Dushort) (int) ((data[p]-b) * one_over_a);
*/
/* MJK 3.24.99 I put this here so if the value is close
   to the missing value and get's rounded up it won't come out
   as missing data */
               if (compdata2[p] == 65535){
                  compdata2[p] = 65534;
               }
            }
         }
         /* TODO: byte-swapping on little endian??? */
#endif
      }
   }

   else {
      /* compressmode==4 */
#ifdef _CRAY
      cray_to_ieee_array( compdata, data, nrncnl );
#else
      /* other machines: just copy 4-byte IEEE floats */
      assert( sizeof(float)==4 );
      memcpy( compdata, data, nrncnl*4 );
      /* TODO: byte-swapping on little endian??? */
#endif
   }
}



/*
 * Decompress a 3-D grid from 1-byte integers to 4-byte floats.
 * Input:  nr, nc, nl - size of grid
 *         compdata - array of [nr*nr*nl*compressmode] bytes
 *         ga, gb - arrays of decompression factors
 *         compressmode - 1, 2 or 4 bytes per grid point
 *         data - address to put decompressed values
 * Output:  data - uncompressed floating point data values
 */
void v5dDecompressGrid( int nr, int nc, int nl, int compressmode,
                        void *compdata, float ga[], float gb[],
                        float data[] )
{
   int nrnc = nr * nc;
   int nrncnl = nr * nc * nl;
   V5Dubyte *compdata1 = (V5Dubyte *) compdata;
   V5Dushort *compdata2 = (V5Dushort *) compdata;

   if (compressmode == 1) {
      int p, i, lev;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float a = ga[lev];
         float b = gb[lev];

         /* WLH 2-2-95 */
         float d, aa;
         int id;
         if (a > 0.0000000001) {
           d = b / a;
           id = floor(d);
           d = d - id;
           aa = a * 0.000001;
         }
         else {
           id = 1;
         }
         if (-254 <= id && id <= 0 && d < aa) {
           for (i=0;i<nrnc;i++,p++) {
              if (compdata1[p]==255) {
                 data[p] = MISSING;
              }
              else {
                 data[p] = (float) (int) compdata1[p] * a + b;
                 if (fabs(data[p]) < aa) data[p] = aa;
              }
           }
         }
         else {
           for (i=0;i<nrnc;i++,p++) {
              if (compdata1[p]==255) {
                 data[p] = MISSING;
              }
              else {
                 data[p] = (float) (int) compdata1[p] * a + b;
              }
           }
         }
         /* end of WLH 2-2-95 */
      }
   }

   else if (compressmode == 2) {
      int p, i, lev;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float a = ga[lev];
         float b = gb[lev];
#ifdef _CRAY
         /* this is tricky because sizeof(V5Dushort)==8, not 2 */
         for (i=0;i<nrnc;i++,p++) {
            int compvalue;
            compvalue = (compdata1[p*2] << 8) | compdata1[p*2+1];
            if (compvalue==65535) {
               data[p] = MISSING;
            }
            else {
               data[p] = (float) compvalue * a + b;
            }
         }
#else
         /* sizeof(V5Dushort)==2! */
         for (i=0;i<nrnc;i++,p++) {
            if (compdata2[p]==65535) {
               data[p] = MISSING;
            }
            else {
               data[p] = (float) (int) compdata2[p] * a + b;
            }
         }
#endif
      }
   }

   else {
      /* compressmode==4 */
#ifdef _CRAY
      ieee_to_cray_array( data, compdata, nrncnl );
#else
      /* other machines: just copy 4-byte IEEE floats */
      assert( sizeof(float)==4 );
      memcpy( data, compdata, nrncnl*4 );
#endif
   }
}




/*
 * Return the size (in bytes) of the 3-D grid specified by time and var.
 * Input:  v - pointer to v5dstruct describing the file
 *         time, var - which timestep and variable
 * Return:  number of data points.
 */
int v5dSizeofGrid( const v5dstruct *v, int time, int var )
{
   return v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
}



/*
 * Initialize a v5dstructure to reasonable initial values.
 * Input:  v - pointer to v5dstruct.
 */
void v5dInitStruct( v5dstruct *v )
{
   int i;

   /* set everything to zero */
   memset( v, 0, sizeof(v5dstruct) );

   /* special cases */
   v->Projection = -1;
   v->VerticalSystem = -1;

   for (i=0;i<MAXVARS;i++) {
      v->MinVal[i] = MISSING;
      v->MaxVal[i] = -MISSING;
      v->LowLev[i] = 0;
   }

   /* set file version */
   strcpy(v->FileVersion, FILE_VERSION);

   v->CompressMode = 1;
   v->FileDesc = -1;
}



/*
 * Return a pointer to a new, initialized v5dstruct.
 */
v5dstruct *v5dNewStruct( void )
{
   v5dstruct *v;

   v = (v5dstruct *) malloc( sizeof(v5dstruct) );
   if (v) {
      v5dInitStruct(v);
   }
   return v;
}



/*
 * Free an initialized v5dstruct. (Todd Plessel)
 */
void v5dFreeStruct( v5dstruct* v )
{
   /*assert( v5dVerifyStruct( v ) );*/
   free( v );
   v = 0;
}



/*
 * Do some checking that the information in a v5dstruct is valid.
 * Input:  v - pointer to v5dstruct
 * Return:  1 = g is ok, 0 = g is invalid
 */
int v5dVerifyStruct( const v5dstruct *v )
{
   int var, i, invalid, maxnl;

   invalid = 0;

   if (!v)
      return 0;

   /* Number of variables */
   if (v->NumVars<0) {
      printf("Invalid number of variables: %d\n", v->NumVars );
      invalid = 1;
   }
   else if (v->NumVars>MAXVARS) {
      printf("Too many variables: %d  (Maximum is %d)\n",
             v->NumVars, MAXVARS);
      invalid = 1;
   }

   /* Variable Names */
   for (i=0;i<v->NumVars;i++) {
      if (v->VarName[i][0]==0) {
         printf("Missing variable name: VarName[%d]=\"\"\n", i );
         invalid = 1;
      }
   }

   /* Number of timesteps */
   if (v->NumTimes<0) {
      printf("Invalid number of timesteps: %d\n", v->NumTimes );
      invalid = 1;
   }
   else if (v->NumTimes>MAXTIMES) {
      printf("Too many timesteps: %d  (Maximum is %d)\n",
             v->NumTimes, MAXTIMES );
      invalid = 1;
   }

   /* Make sure timestamps are increasing */
   for (i=1;i<v->NumTimes;i++) {
      int date0 = v5dYYDDDtoDays( v->DateStamp[i-1] );
      int date1 = v5dYYDDDtoDays( v->DateStamp[i] );
      int time0 = v5dHHMMSStoSeconds( v->TimeStamp[i-1] );
      int time1 = v5dHHMMSStoSeconds( v->TimeStamp[i] );
      if (time1<=time0 && date1<=date0) {
         printf("Timestamp for step %d must be later than step %d\n", i, i-1);
         invalid = 1;
      }
   }

   /* Rows */
   if (v->Nr<2) {
      printf("Too few rows: %d (2 is minimum)\n", v->Nr );
      invalid = 1;
   }
   else if (v->Nr>MAXROWS) {
      printf("Too many rows: %d (%d is maximum)\n", v->Nr, MAXROWS );
      invalid = 1;
   }

   /* Columns */
   if (v->Nc<2) {
      printf("Too few columns: %d (2 is minimum)\n", v->Nc );
      invalid = 1;
   }
   else if (v->Nc>MAXCOLUMNS) {
      printf("Too many columns: %d (%d is maximum)\n", v->Nc, MAXCOLUMNS );
      invalid = 1;
   }

   /* Levels */
   maxnl = 0;
   for (var=0;var<v->NumVars;var++) {
      if (v->LowLev[var] < 0) {
         printf("Low level cannot be negative for var %s: %d\n",
                 v->VarName[var], v->LowLev[var] );
         invalid = 1;
      }
      if (v->Nl[var]<1) {
         printf("Too few levels for var %s: %d (1 is minimum)\n",
                 v->VarName[var], v->Nl[var] );
         invalid = 1;
      }
      if (v->Nl[var]+v->LowLev[var]>MAXLEVELS) {
         printf("Too many levels for var %s: %d (%d is maximum)\n",
                 v->VarName[var], v->Nl[var]+v->LowLev[var], MAXLEVELS );
         invalid = 1;
      }
      if (v->Nl[var]+v->LowLev[var]>maxnl) {
         maxnl = v->Nl[var]+v->LowLev[var];
      }
   }

   if (v->CompressMode != 1 && v->CompressMode != 2 && v->CompressMode != 4) {
      printf("Bad CompressMode: %d (must be 1, 2 or 4)\n", v->CompressMode );
      invalid = 1;
   }

   switch (v->VerticalSystem) {
      case 0:
      case 1:
         if (v->VertArgs[1]==0.0) {
            printf("Vertical level increment is zero, must be non-zero\n");
            invalid = 1;
         }
         break;
      case 2:
         /* Check that Height values increase upward */
         for (i=1;i<maxnl;i++) {
            if (v->VertArgs[i] <= v->VertArgs[i-1]) {
               printf("Height[%d]=%f <= Height[%d]=%f, level heights must increase\n",
                      i, v->VertArgs[i], i-1, v->VertArgs[i-1] );
               invalid = 1;
               break;
            }
         }
         break;
      case 3:
         /* Check that Pressure values decrease upward */
         for (i=1;i<maxnl;i++) {
            if (v->VertArgs[i] <= v->VertArgs[i-1]) {
               printf("Pressure[%d]=%f >= Pressure[%d]=%f, level pressures must decrease\n",
                      i, height_to_pressure(v->VertArgs[i]),
                      i-1, height_to_pressure(v->VertArgs[i-1]) );
               invalid = 1;
               break;
            }
         }
         break;
      default:
         printf("VerticalSystem = %d, must be in 0..3\n", v->VerticalSystem );
         invalid = 1;
   }


   switch (v->Projection) {
      case 0:  /* Generic */
         if (v->ProjArgs[2]==0.0) {
            printf("Row Increment (ProjArgs[2]) can't be zero\n");
            invalid = 1;
         }
         if (v->ProjArgs[3]==0.0) {
            printf("Column increment (ProjArgs[3]) can't be zero\n");
            invalid = 1;
         }
         break;
      case 1:  /* Cylindrical equidistant */
         if (v->ProjArgs[2]<0.0) {
            printf("Row Increment (ProjArgs[2]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[2] );
            invalid = 1;
         }
         if (v->ProjArgs[3]<=0.0) {
            printf("Column Increment (ProjArgs[3]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[3] );
            invalid = 1;
         }
         break;
      case 2:  /* Lambert Conformal */
         if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
            printf("Lat1 (ProjArgs[0]) out of range: %g\n", v->ProjArgs[0] );
            invalid = 1;
         }
         if (v->ProjArgs[1]<-90.0 || v->ProjArgs[1]>90.0) {
            printf("Lat2 (ProjArgs[1] out of range: %g\n", v->ProjArgs[1] );
            invalid = 1;
         }
         if (v->ProjArgs[5]<=0.0) {
            printf("ColInc (ProjArgs[5]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[5] );
            invalid = 1;
         }
         break;
      case 3:  /* Stereographic */
         if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
            printf("Central Latitude (ProjArgs[0]) out of range: ");
            printf("%g  (must be in +/-90)\n", v->ProjArgs[0] );
            invalid = 1;
         }
         if (v->ProjArgs[1]<-180.0 || v->ProjArgs[1]>180.0) {
            printf("Central Longitude (ProjArgs[1]) out of range: ");
            printf("%g  (must be in +/-180)\n", v->ProjArgs[1] );
            invalid = 1;
         }
         if (v->ProjArgs[4]<0) {
            printf("Column spacing (ProjArgs[4]) = %g  (must be positive)\n",
                   v->ProjArgs[4]);
            invalid = 1;
         }
         break;
      case 4:  /* Rotated */
         /* WLH 4-21-95 */
         if (v->ProjArgs[2]<=0.0) {
            printf("Row Increment (ProjArgs[2]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[2] );
            invalid = 1;
         }
         if (v->ProjArgs[3]<=0.0) {
            printf("Column Increment = (ProjArgs[3]) %g  (must be >=0.0)\n",
                   v->ProjArgs[3] );
            invalid = 1;
         }
         if (v->ProjArgs[4]<-90.0 || v->ProjArgs[4]>90.0) {
            printf("Central Latitude (ProjArgs[4]) out of range: ");
            printf("%g  (must be in +/-90)\n", v->ProjArgs[4] );
            invalid = 1;
         }
         if (v->ProjArgs[5]<-180.0 || v->ProjArgs[5]>180.0) {
            printf("Central Longitude (ProjArgs[5]) out of range: ");
            printf("%g  (must be in +/-180)\n", v->ProjArgs[5] );
            invalid = 1;
         }
         if (v->ProjArgs[6]<-180.0 || v->ProjArgs[6]>180.0) {
            printf("Central Longitude (ProjArgs[6]) out of range: ");
            printf("%g  (must be in +/-180)\n", v->ProjArgs[6] );
            invalid = 1;
         }
         break;
      case 5: /*Mercator*/
         if (v->ProjArgs[2] == 0.0){
            printf("Row Increment(Km) can not be 0.0\n");
            invalid = 1;
         }
         if (v->ProjArgs[3] == 0.0){
            printf("Column Increment(Km) can not be 0.0\n");
            invalid = 1;
         }
         break;
      default:
         printf("Projection = %d, must be in 0..4\n", v->Projection );
         invalid = 1;
   }

   return !invalid;
}



/*
 * Get the McIDAS file number and grid number associated with the grid
 * identified by time and var.
 * Input:  v - v5d grid struct
 *         time, var - timestep and variable of grid
 * Output:  mcfile, mcgrid - McIDAS grid file number and grid number
 */
int v5dGetMcIDASgrid( v5dstruct *v, int time, int var,
                      int *mcfile, int *mcgrid )
{
   if (time<0 || time>=v->NumTimes) {
      printf("Bad time argument to v5dGetMcIDASgrid: %d\n", time );
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Bad var argument to v5dGetMcIDASgrid: %d\n", var );
      return 0;
   }

   *mcfile = (int) v->McFile[time][var];
   *mcgrid = (int) v->McGrid[time][var];
   return 1;
}



/*
 * Set the McIDAS file number and grid number associated with the grid
 * identified by time and var.
 * Input:  v - v5d grid struct
 *         time, var - timestep and variable of grid
 *         mcfile, mcgrid - McIDAS grid file number and grid number
 * Return:  1 = ok, 0 = error (bad time or var)
 */
int v5dSetMcIDASgrid( v5dstruct *v, int time, int var,
                      int mcfile, int mcgrid )
{
   if (time<0 || time>=v->NumTimes) {
      printf("Bad time argument to v5dSetMcIDASgrid: %d\n", time );
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Bad var argument to v5dSetMcIDASgrid: %d\n", var );
      return 0;
   }

   v->McFile[time][var] = (short) mcfile;
   v->McGrid[time][var] = (short) mcgrid;
   return 1;
}



/**********************************************************************/
/*****                    Input Functions                         *****/
/**********************************************************************/



/*
 * Read the header from a COMP* file and return results in the v5dstruct.
 * Input:  f - the file descriptor
 *         v - pointer to a v5dstruct.
 * Return:  1 = ok, 0 = error.
 */
static int read_comp_header( int f, v5dstruct *v )
{
   unsigned int id;

   /* reset file position to start of file */
   lseek( f, 0, SEEK_SET );

   /* read file ID */
   read_int4( f, (int *) &id );

   if (id==0x80808080 || id==0x80808081) {
      /* Older COMP5D format */
      int gridtimes, gridparms;
      int i, j, it, iv, nl;
      int gridsize;
      float hgttop, hgtinc;
      /*char *compgrid;*/

      if (id==0x80808080) {
         /* 20 vars, 300 times */
         gridtimes = 300;
         gridparms = 20;
      }
      else {
         /* 30 vars, 400 times */
         gridtimes = 400;
         gridparms = 30;
      }

      v->FirstGridPos = 12*4 + 8*gridtimes + 4*gridparms;

      read_int4( f, &v->NumTimes );
      read_int4( f, &v->NumVars );
      read_int4( f, &v->Nr );
      read_int4( f, &v->Nc );
      read_int4( f, &nl );
      for (i=0;i<v->NumVars;i++) {
         v->Nl[i] = nl;
         v->LowLev[i] = 0;
      }
      read_float4( f, &v->ProjArgs[0] );
      read_float4( f, &v->ProjArgs[1] );
      read_float4( f, &hgttop );
      read_float4( f, &v->ProjArgs[2] );
      read_float4( f, &v->ProjArgs[3] );
      read_float4( f, &hgtinc );
/*
      for (i=0;i<nl;i++) {
         v->Height[nl-i-1] = hgttop - i * hgtinc;
      }
*/
      v->VerticalSystem = 1;
      v->VertArgs[0] = hgttop - hgtinc * (nl-1);
      v->VertArgs[1] = hgtinc;

      /* read dates and times */
      for (i=0;i<gridtimes;i++) {
         read_int4( f, &j );
         v->DateStamp[i] = v5dDaysToYYDDD( j );
      }
      for (i=0;i<gridtimes;i++) {
         read_int4( f, &j );
         v->TimeStamp[i] = v5dSecondsToHHMMSS( j );
      }

      /* read variable names */
      for (i=0;i<gridparms;i++) {
         char name[4];
         read_bytes( f, name, 4 );
         /* remove trailing spaces, if any */
         for (j=3;j>0;j--) {
            if (name[j]==' ' || name[j]==0)
              name[j] = 0;
            else
              break;
         }
         strncpy( v->VarName[i], name, 4 );
         v->VarName[i][4] = 0;
      }

      gridsize = ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
      for (i=0;i<v->NumVars;i++) {
         v->GridSize[i] = 8 + gridsize;
      }
      v->SumGridSizes = (8+gridsize) * v->NumVars;

      /* read the grids and their ga,gb values to find min and max values */

      for (i=0;i<v->NumVars;i++) {
         v->MinVal[i] = 999999.9;
         v->MaxVal[i] = -999999.9;
      }

      /*compgrid = (char *) malloc( gridsize );*/

      for (it=0; it<v->NumTimes; it++) {
         for (iv=0; iv<v->NumVars; iv++) {
            float ga, gb;
            float min, max;

            read_float4( f, &ga );
            read_float4( f, &gb );

            /* skip ahead by 'gridsize' bytes */
            if (lseek( f, gridsize, SEEK_CUR )==-1) {
               printf("Error:  Unexpected end of file, ");
               printf("file may be corrupted.\n");
               return 0;
            }
            min = -(125.0+gb)/ga;
            max = (125.0-gb)/ga;
            if (min<v->MinVal[iv])  v->MinVal[iv] = min;
            if (max>v->MaxVal[iv])  v->MaxVal[iv] = max;
         }
      }

      /*free( compgrid );*/

      /* done */
   }
   else if (id==0x80808082 || id==0x80808083) {
      /* Newer COMP5D format */
      int gridtimes, gridsize;
      int it, iv, nl, i, j;
      float delta;

      read_int4( f, &gridtimes );
      read_int4( f, &v->NumVars );
      read_int4( f, &v->NumTimes );
      read_int4( f, &v->Nr );
      read_int4( f, &v->Nc );
      read_int4( f, &nl );
      for (i=0;i<v->NumVars;i++) {
         v->Nl[i] = nl;
      }

      read_float4( f, &v->ProjArgs[2] );
      read_float4( f, &v->ProjArgs[3] );

      /* Read height and determine if equal spacing */
      v->VerticalSystem = 1;
      for (i=0;i<nl;i++) {
         read_float4( f, &v->VertArgs[i] );
         if (i==1) {
            delta = v->VertArgs[1] - v->VertArgs[0];
         }
         else if (i>1) {
            if (delta != (v->VertArgs[i] - v->VertArgs[i-1])) {
               v->VerticalSystem = 2;
            }
         }
      }
      if (v->VerticalSystem==1) {
         v->VertArgs[1] = delta;
      }

      /* read variable names */
      for (iv=0; iv<v->NumVars; iv++) {
         char name[8];

         read_bytes( f, name, 8 );

         /* remove trailing spaces, if any */
         for (j=7;j>0;j--) {
            if (name[j]==' ' || name[j]==0)
              name[j] = 0;
            else
              break;
         }
         strncpy( v->VarName[iv], name, 8 );
         v->VarName[iv][8] = 0;
      }

      for (iv=0;iv<v->NumVars;iv++) {
         read_float4( f, &v->MinVal[iv] );
      }
      for (iv=0;iv<v->NumVars;iv++) {
         read_float4( f, &v->MaxVal[iv] );
      }
      for (it=0;it<gridtimes;it++) {
         read_int4( f, &j );
         v->TimeStamp[it] = v5dSecondsToHHMMSS( j );
      }
      for (it=0;it<gridtimes;it++) {
         read_int4( f, &j );
         v->DateStamp[it] = v5dDaysToYYDDD( j );
      }
      for (it=0;it<gridtimes;it++) {
         float nlat;
         read_float4( f, &nlat );
         if (it==0)  v->ProjArgs[0] = nlat;
      }
      for (it=0;it<gridtimes;it++) {
         float wlon;
         read_float4( f, &wlon );
         if (it==0)  v->ProjArgs[1] = wlon;
      }

      /* calculate grid storage sizes */
      if (id==0x80808082) {
         gridsize = nl*2*4 + ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
      }
      else {
         /* McIDAS grid and file numbers present */
         gridsize = 8 + nl*2*4 + ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
      }
      for (i=0;i<v->NumVars;i++) {
         v->GridSize[i] = gridsize;
      }
      v->SumGridSizes = gridsize * v->NumVars;

      /* read McIDAS numbers??? */

      /* size (in bytes) of all header info */
      v->FirstGridPos = 9*4 + v->Nl[0]*4 + v->NumVars*16 + gridtimes*16;

   }

   v->CompressMode = 1; /* one byte per grid point */
   v->Projection = 1;  /* Cylindrical equidistant */
   v->FileVersion[0] = 0;

   return 1;
}



/*
 * Read a compressed grid from a COMP* file.
 * Return:  1 = ok, 0 = error.
 */
static int read_comp_grid( v5dstruct *v, int time, int var,
                           float *ga, float *gb, void *compdata )
{
   unsigned int pos;
   V5Dubyte bias;
   int i, n, nl;
   int f;
   V5Dubyte *compdata1 = (V5Dubyte *) compdata;

   f = v->FileDesc;

   /* move to position in file */
   pos = grid_position( v, time, var );
   lseek( f, pos, SEEK_SET );

   if (v->FileFormat==0x80808083) {
      /* read McIDAS grid and file numbers */
      int mcfile, mcgrid;
      read_int4( f, &mcfile );
      read_int4( f, &mcgrid );
      v->McFile[time][var] = (short) mcfile;
      v->McGrid[time][var] = (short) mcgrid;
   }

   nl = v->Nl[var];

   if (v->FileFormat==0x80808080 || v->FileFormat==0x80808081) {
      /* single ga,gb pair for whole grid */
      float a, b;
      read_float4( f, &a );
      read_float4( f, &b );
      /* convert a, b to new v5d ga, gb values */
      for (i=0;i<nl;i++) {
         if (a==0.0) {
            ga[i] = gb[i] = 0.0;
         }
         else {
            gb[i] = (b+128.0) / -a;
            ga[i] = 1.0 / a;
         }
      }
      bias = 128;
   }
   else {
      /* read ga, gb arrays */
      read_float4_array( f, ga, v->Nl[var] );
      read_float4_array( f, gb, v->Nl[var] );

      /* convert ga, gb values to v5d system */
      for (i=0;i<nl;i++) {
         if (ga[i]==0.0) {
            ga[i] = gb[i] = 0.0;
         }
         else {
            /*gb[i] = (gb[i]+125.0) / -ga[i];*/
            gb[i] = (gb[i]+128.0) / -ga[i];
            ga[i] = 1.0 / ga[i];
         }
      }
      bias = 128;  /* 125 ??? */
   }

   /* read compressed grid data */
   n = v->Nr * v->Nc * v->Nl[var];
   if (read_bytes( f, compdata1, n )!=n)
      return 0;

   /* convert data values to v5d system */
   n = v->Nr * v->Nc * v->Nl[var];
   for (i=0;i<n;i++) {
      compdata1[i] += bias;
   }

   return 1;
}



/*
 * Read a v5d file header.
 * Input:  f - file opened for reading.
 *         v - pointer to v5dstruct to store header info into.
 * Return:  1 = ok, 0 = error.
 */
static int read_v5d_header( v5dstruct *v )
{
#define SKIP(N)   lseek( f, N, SEEK_CUR )
   int end_of_header = 0;
   unsigned int id;
   int idlen, var, numargs;
   int f;

   f = v->FileDesc;

   /* first try to read the header id */
   read_int4( f, (int*) &id );
   read_int4( f, &idlen );
   if (id==TAG_ID && idlen==0) {
      /* this is a v5d file */
      v->FileFormat = 0;
   }
   else if (id>=0x80808080 && id<=0x80808083) {
      /* this is an old COMP* file */
      v->FileFormat = id;
      return read_comp_header( f, v );
   }
   else {
      /* unknown file type */
      printf("Error: not a v5d file\n");
      return 0;
   }

   v->CompressMode = 1; /* default */

   while (!end_of_header) {
      int tag, length;
      int i, var, time, nl, lev;

      if (read_int4(f,&tag)<1 || read_int4(f,&length)<1) {
         printf("Error while reading header, premature EOF\n");
         return 0;
      }

      switch (tag) {
         case TAG_VERSION:
            assert( length==10 );
            read_bytes( f, v->FileVersion, 10 );
            /* Check if reading a file made by a future version of Vis5D */
            if (strcmp(v->FileVersion, FILE_VERSION)>0) {
               /* WLH 6 Oct 98 */
               printf("Warning: Trying to read a version %s file,", v->FileVersion);
               printf(" you should upgrade Vis5D.\n");
            }
            break;
         case TAG_NUMTIMES:
            assert( length==4 );
            read_int4( f, &v->NumTimes );
            break;
         case TAG_NUMVARS:
            assert( length==4 );
            read_int4( f, &v->NumVars );
            break;
         case TAG_VARNAME:
            assert( length==14 );   /* 1 int + 10 char */
            read_int4( f, &var );
            read_bytes( f, v->VarName[var], 10 );
            break;
         case TAG_NR:
            /* Number of rows for all variables */
            assert( length==4 );
            read_int4( f, &v->Nr );
            break;
         case TAG_NC:
            /* Number of columns for all variables */
            assert( length==4 );
            read_int4( f, &v->Nc );
            break;
         case TAG_NL:
            /* Number of levels for all variables */
            assert( length==4 );
            read_int4( f, &nl );
            for (i=0;i<v->NumVars;i++) {
               v->Nl[i] = nl;
            }
            break;
         case TAG_NL_VAR:
            /* Number of levels for one variable */
            assert( length==8 );
            read_int4( f, &var );
            read_int4( f, &v->Nl[var] );
            break;
         case TAG_LOWLEV_VAR:
            /* Lowest level for one variable */
            assert( length==8 );
            read_int4( f, &var );
            read_int4( f, &v->LowLev[var] );
            break;

         case TAG_TIME:
            /* Time stamp for 1 timestep */
            assert( length==8 );
            read_int4( f, &time );
            read_int4( f, &v->TimeStamp[time] );
            break;
         case TAG_DATE:
            /* Date stamp for 1 timestep */
            assert( length==8 );
            read_int4( f, &time );
            read_int4( f, &v->DateStamp[time] );
            break;

         case TAG_MINVAL:
            /* Minimum value for a variable */
            assert( length==8 );
            read_int4( f, &var );
            read_float4( f, &v->MinVal[var] );
            break;
         case TAG_MAXVAL:
            /* Maximum value for a variable */
            assert( length==8 );
            read_int4( f, &var );
            read_float4( f, &v->MaxVal[var] );
            break;
         case TAG_COMPRESS:
            /* Compress mode */
            assert( length==4 );
            read_int4( f, &v->CompressMode );
            break;
         case TAG_UNITS:
            /* physical units */
            assert( length==24 );
            read_int4( f, &var );
            read_bytes( f, v->Units[var], 20 );
            break;

         /*
          * Vertical coordinate system
          */
         case TAG_VERTICAL_SYSTEM:
            assert( length==4 );
            read_int4( f, &v->VerticalSystem );
            if (v->VerticalSystem<0 || v->VerticalSystem>3) {
               printf("Error: bad vertical coordinate system: %d\n",
                      v->VerticalSystem );
            }
            break;
         case TAG_VERT_ARGS:
            read_int4( f, &numargs );
            assert( numargs <= MAXVERTARGS );
            read_float4_array( f, v->VertArgs, numargs );
            assert( length==numargs*4+4 );
            break;
         case TAG_HEIGHT:
            /* height of a grid level */
            assert( length==8 );
            read_int4( f, &lev );
            read_float4( f, &v->VertArgs[lev] );
            break;
         case TAG_BOTTOMBOUND:
            assert( length==4 );
            read_float4( f, &v->VertArgs[0] );
            break;
         case TAG_LEVINC:
            assert( length==4 );
            read_float4( f, &v->VertArgs[1] );
            break;

         /*
          * Map projection information
          */
         case TAG_PROJECTION:
            assert( length==4 );
            read_int4( f, &v->Projection );
            if (v->Projection<0 || v->Projection>5) { /* WLH 4-21-95 */
               printf("Error while reading header, bad projection (%d)\n",
                       v->Projection );
               return 0;
            }
            break;
         case TAG_PROJ_ARGS:
            read_int4( f, &numargs );
            assert( numargs <= MAXPROJARGS );
            read_float4_array( f, v->ProjArgs, numargs );
            assert( length==4*numargs+4 );
            break;
         case TAG_NORTHBOUND:
            assert( length==4 );
            if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
               read_float4( f, &v->ProjArgs[0] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_WESTBOUND:
            assert( length==4 );
            if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
               read_float4( f, &v->ProjArgs[1] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_ROWINC:
            assert( length==4 );
            if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
               read_float4( f, &v->ProjArgs[2] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_COLINC:
            assert( length==4 );
            if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
               read_float4( f, &v->ProjArgs[3] );
            }
            else if (v->Projection==2) {
               read_float4( f, &v->ProjArgs[5] );
            }
            else if (v->Projection==3) {
               read_float4( f, &v->ProjArgs[4] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_ROWINCKM:
            assert( length==4 );
            if (v->Projection==5){
               read_float4( f, &v->ProjArgs[2] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_COLINCKM:
            assert( length==4 );
            if (v->Projection==5){
               read_float4( f, &v->ProjArgs[3] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_LAT1:
            assert( length==4 );
            if (v->Projection==2) {
               read_float4( f, &v->ProjArgs[0] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_LAT2:
            assert( length==4 );
            if (v->Projection==2) {
               read_float4( f, &v->ProjArgs[1] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_POLE_ROW:
            assert( length==4 );
            if (v->Projection==2) {
               read_float4( f, &v->ProjArgs[2] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_POLE_COL:
            assert( length==4 );
            if (v->Projection==2) {
               read_float4( f, &v->ProjArgs[3] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_CENTLON:
            assert( length==4 );
            if (v->Projection==2) {
               read_float4( f, &v->ProjArgs[4] );
            }
            else if (v->Projection==3 || v->Projection==5) {
               read_float4( f, &v->ProjArgs[1] );
            }
            else if (v->Projection==4) { /* WLH 4-21-95 */
               read_float4( f, &v->ProjArgs[5] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_CENTLAT:
            assert( length==4 );
            if (v->Projection==3 || v->Projection==5) {
               read_float4( f, &v->ProjArgs[0] );
            }
            else if (v->Projection==4) { /* WLH 4-21-95 */
               read_float4( f, &v->ProjArgs[4] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_CENTROW:
            assert( length==4 );
            if (v->Projection==3) {
               read_float4( f, &v->ProjArgs[2] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_CENTCOL:
            assert( length==4 );
            if (v->Projection==3) {
               read_float4( f, &v->ProjArgs[3] );
            }
            else {
               SKIP( 4 );
            }
            break;
         case TAG_ROTATION:
            assert( length==4 );
            if (v->Projection==4) { /* WLH 4-21-95 */
               read_float4( f, &v->ProjArgs[6] );
            }
            else {
               SKIP( 4 );
            }
            break;

         case TAG_END:
            /* end of header */
            end_of_header = 1;
            lseek( f, length, SEEK_CUR );
            break;

         default:
            /* unknown tag, skip to next tag */
            printf("Unknown tag: %d  length=%d\n", tag, length );
            lseek( f, length, SEEK_CUR );
            break;
      }

   }

   v5dVerifyStruct( v );

   /* Now we're ready to read the grid data */

   /* Save current file pointer */
   v->FirstGridPos = ltell(f);

   /* compute grid sizes */
   v->SumGridSizes = 0;
   for (var=0;var<v->NumVars;var++) {
      v->GridSize[var] = 8 * v->Nl[var] + v5dSizeofGrid( v, 0, var );
      v->SumGridSizes += v->GridSize[var];
   }

   return 1;
#undef SKIP
}




/*
 * Open a v5d file for reading.
 * Input:  filename - name of v5d file to open
 *         v - pointer to a v5dstruct in which to put header info or NULL
 *             if a struct should be dynamically allocated.
 * Return:  NULL if error, else v or a pointer to a new v5dstruct if v was NULL
 */
v5dstruct *v5dOpenFile( const char *filename, v5dstruct *v )
{
   int fd;

   fd = open( filename, O_RDONLY );
   if (fd==-1) {
      /* error */
      return 0;
   }

   if (v) {
      v5dInitStruct( v );
   }
   else {
      v = v5dNewStruct();
      if (!v) {
         return NULL;
      }
   }

   v->FileDesc = fd;
   v->Mode = 'r';
   if (read_v5d_header( v )) {
      return v;
   }
   else {
      return NULL;
   }
}




/*
 * Read a compressed grid from a v5d file.
 * Input:  v - pointer to v5dstruct describing the file
 *         time, var - which timestep and variable
 *         ga, gb - arrays to store grid (de)compression values
 *         compdata - address of where to store compressed grid data.
 * Return:  1 = ok, 0 = error.
 */
int v5dReadCompressedGrid( v5dstruct *v, int time, int var,
                           float *ga, float *gb, void *compdata )
{
   int pos, n, k;

   if (time<0 || time>=v->NumTimes) {
      printf("Error in v5dReadCompressedGrid: bad timestep argument (%d)\n",
             time);
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Error in v5dReadCompressedGrid: bad var argument (%d)\n",
             var);
      return 0;
   }

   if (v->FileFormat) {
      /* old COMP* file */
      return read_comp_grid( v, time, var, ga, gb, compdata );
   }

   /* move to position in file */
   pos = grid_position( v, time, var );
   lseek( v->FileDesc, pos, SEEK_SET );

   /* read ga, gb arrays */
   read_float4_array( v->FileDesc, ga, v->Nl[var] );
   read_float4_array( v->FileDesc, gb, v->Nl[var] );

   /* read compressed grid data */
   n = v->Nr * v->Nc * v->Nl[var];
   if (v->CompressMode==1) {
      k = read_block( v->FileDesc, compdata, n, 1 )==n;
   }
   else if (v->CompressMode==2) {
      k = read_block( v->FileDesc, compdata, n, 2 )==n;
   }
   else if (v->CompressMode==4) {
      k = read_block( v->FileDesc, compdata, n, 4 )==n;
   }
   if (!k) {
      /* error */
      printf("Error in v5dReadCompressedGrid: read failed, bad file?\n");
   }
   return k;


/*
   n = v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
   if (read( v->FileDesc, compdata, n )==n)
      return 1;
   else
      return 0;
*/
}




/*
 * Read a grid from a v5d file, decompress it and return it.
 * Input:  v - pointer to v5dstruct describing file header
 *         time, var - which timestep and variable.
 *         data - address of buffer to put grid data
 * Output:  data - the grid data
 * Return:  1 = ok, 0 = error.
 */
int v5dReadGrid( v5dstruct *v, int time, int var, float data[] )
{
   float ga[MAXLEVELS], gb[MAXLEVELS];
   void *compdata;
   int bytes;

   if (time<0 || time>=v->NumTimes) {
      printf("Error in v5dReadGrid: bad timestep argument (%d)\n", time);
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Error in v5dReadGrid: bad variable argument (%d)\n", var);
      return 0;
   }

   /* allocate compdata buffer */
   if (v->CompressMode==1) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned char);
   }
   else if (v->CompressMode==2) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned short);
   }
   else if (v->CompressMode==4) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(float);
   }
   compdata = (void *) malloc( bytes );
   if (!compdata) {
      printf("Error in v5dReadGrid: out of memory (needed %d bytes)\n", bytes);
      return 0;
   }

   /* read the compressed data */
   if (!v5dReadCompressedGrid( v, time, var, ga, gb, compdata )) {
      return 0;
   }

   /* decompress the data */
   v5dDecompressGrid( v->Nr, v->Nc, v->Nl[var], v->CompressMode,
                      compdata, ga, gb, data );

   /* free compdata */
   free( compdata );
   return 1;
}




/**********************************************************************/
/*****                   Output Functions                         *****/
/**********************************************************************/



static int write_tag( v5dstruct *v, int tag, int length, int newfile )
{
   if (!newfile) {
      /* have to check that there's room in header to write this tagged item */
      if (v->CurPos+8+length > v->FirstGridPos) {
         printf("Error: out of header space!\n");
         /* Out of header space! */
         return 0;
      }
   }

   if (write_int4( v->FileDesc, tag )==0)  return 0;
   if (write_int4( v->FileDesc, length )==0)  return 0;
   v->CurPos += 8 + length;
   return 1;
}



/*
 * Write the information in the given v5dstruct as a v5d file header.
 * Note that the current file position is restored when this function
 * returns normally.
 * Input:  f - file already open for writing
 *         v - pointer to v5dstruct
 * Return:  1 = ok, 0 = error.
 */
static int write_v5d_header( v5dstruct *v )
{
   int var, time, filler, maxnl;
   int f;
   int newfile;

   if (v->FileFormat!=0) {
      printf("Error: v5d library can't write comp5d format files.\n");
      return 0;
   }

   f = v->FileDesc;

   if (!v5dVerifyStruct( v ))
      return 0;

   /* Determine if we're writing to a new file */
   if (v->FirstGridPos==0) {
      newfile = 1;
   }
   else {
      newfile = 0;
   }

   /* compute grid sizes */
   v->SumGridSizes = 0;
   for (var=0;var<v->NumVars;var++) {
      v->GridSize[var] = 8 * v->Nl[var] + v5dSizeofGrid( v, 0, var );
      v->SumGridSizes += v->GridSize[var];
   }

   /* set file pointer to start of file */
   lseek( f, 0, SEEK_SET );
   v->CurPos = 0;

   /*
    * Write the tagged header info
    */
#define WRITE_TAG( V, T, L )  if (!write_tag(V,T,L,newfile))  return 0;

   /* ID */
   WRITE_TAG( v, TAG_ID, 0 );

   /* File Version */
   WRITE_TAG( v, TAG_VERSION, 10 );
   write_bytes( f, FILE_VERSION, 10 );

   /* Number of timesteps */
   WRITE_TAG( v, TAG_NUMTIMES, 4 );
   write_int4( f, v->NumTimes );

   /* Number of variables */
   WRITE_TAG( v, TAG_NUMVARS, 4 );
   write_int4( f, v->NumVars );

   /* Names of variables */
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_VARNAME, 14 );
      write_int4( f, var );
      write_bytes( f, v->VarName[var], 10 );
   }

   /* Physical Units */
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_UNITS, 24 );
      write_int4( f, var );
      write_bytes( f, v->Units[var], 20 );
   }

   /* Date and time of each timestep */
   for (time=0;time<v->NumTimes;time++) {
      WRITE_TAG( v, TAG_TIME, 8 );
      write_int4( f, time );
      write_int4( f, v->TimeStamp[time] );
      WRITE_TAG( v, TAG_DATE, 8 );
      write_int4( f, time );
      write_int4( f, v->DateStamp[time] );
   }

   /* Number of rows */
   WRITE_TAG( v, TAG_NR, 4 );
   write_int4( f, v->Nr );

   /* Number of columns */
   WRITE_TAG( v, TAG_NC, 4 );
   write_int4( f, v->Nc );

   /* Number of levels, compute maxnl */
   maxnl = 0;
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_NL_VAR, 8 );
      write_int4( f, var );
      write_int4( f, v->Nl[var] );
      WRITE_TAG( v, TAG_LOWLEV_VAR, 8 );
      write_int4( f, var );
      write_int4( f, v->LowLev[var] );
      if (v->Nl[var]+v->LowLev[var]>maxnl) {
         maxnl = v->Nl[var]+v->LowLev[var];
      }
   }

   /* Min/Max values */
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_MINVAL, 8 );
      write_int4( f, var );
      write_float4( f, v->MinVal[var] );
      WRITE_TAG( v, TAG_MAXVAL, 8 );
      write_int4( f, var );
      write_float4( f, v->MaxVal[var] );
   }

   /* Compress mode */
   WRITE_TAG( v, TAG_COMPRESS, 4 );
   write_int4( f, v->CompressMode );

   /* Vertical Coordinate System */
   WRITE_TAG( v, TAG_VERTICAL_SYSTEM, 4 );
   write_int4( f, v->VerticalSystem );
   WRITE_TAG( v, TAG_VERT_ARGS, 4+4*MAXVERTARGS );
   write_int4( f, MAXVERTARGS );
   write_float4_array( f, v->VertArgs, MAXVERTARGS );

   /* Map Projection */
   WRITE_TAG( v, TAG_PROJECTION, 4 );
   write_int4( f, v->Projection );
   WRITE_TAG( v, TAG_PROJ_ARGS, 4+4*MAXPROJARGS );
   write_int4( f, MAXPROJARGS );
   write_float4_array( f, v->ProjArgs, MAXPROJARGS );

   /* write END tag */
   if (newfile) {
      /* We're writing to a brand new file.  Reserve 10000 bytes */
      /* for future header growth. */
      WRITE_TAG( v, TAG_END, 10000 );
      lseek( f, 10000, SEEK_CUR );

      /* Let file pointer indicate where first grid is stored */
      v->FirstGridPos = ltell( f );
   }
   else {
      /* we're rewriting a header */
      filler = v->FirstGridPos - ltell(f);
      WRITE_TAG( v, TAG_END, filler-8 );
   }

#undef WRITE_TAG

   return 1;
}



/*
 * Open a v5d file for writing.  If the named file already exists,
 * it will be deleted.
 * Input:  filename - name of v5d file to create.
 *         v - pointer to v5dstruct with the header info to write.
 * Return:  1 = ok, 0 = error.
 */
int v5dCreateFile( const char *filename, v5dstruct *v )
{
   mode_t mask;
   int fd;

   mask = 0666;
   fd = open( filename, O_WRONLY | O_CREAT | O_TRUNC, mask );
   if (fd==-1) {
      printf("Error in v5dCreateFile: open failed\n");
      v->FileDesc = -1;
      v->Mode = 0;
      return 0;
   }
   else {
      /* ok */
      v->FileDesc = fd;
      v->Mode = 'w';
      /* write header and return status */
      return write_v5d_header(v);
   }
}



/*
 * Open a v5d file for updating/appending and read the header info.
 * Input:  filename - name of v5d file to open for updating.
 *         v - pointer to v5dstruct in which the file header info will be
 *             put.  If v is NULL a v5dstruct will be allocated and returned.
 * Return:  NULL if error, else v or a pointer to a new v5dstruct if v as NULL
 */
v5dstruct *v5dUpdateFile( const char *filename, v5dstruct *v )
{
   int fd;

   fd = open( filename, O_RDWR );
   if (fd==-1) {
      return NULL;
   }

   if (!v) {
      v = v5dNewStruct();
      if (!v) {
         return NULL;
      }
   }

   v->FileDesc = fd;
   v->Mode = 'w';

   if (read_v5d_header( v )) {
      return v;
   }
   else {
      return NULL;
   }
}



/*
 * Write a compressed grid to a v5d file.
 * Input:  v - pointer to v5dstruct describing the file
 *         time, var - which timestep and variable
 *         ga, gb - the GA and GB (de)compression value arrays
 *         compdata - address of array of compressed data values
 * Return:  1 = ok, 0 = error.
 */
int v5dWriteCompressedGrid( const v5dstruct *v, int time, int var,
                            const float *ga, const float *gb,
                            const void *compdata )
{
   int pos, n, k;

   /* simple error checks */
   if (v->Mode!='w') {
      printf("Error in v5dWriteCompressedGrid: file opened for reading,");
      printf(" not writing.\n");
      return 0;
   }
   if (time<0 || time>=v->NumTimes) {
      printf("Error in v5dWriteCompressedGrid: bad timestep argument (%d)\n",
             time);
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Error in v5dWriteCompressedGrid: bad variable argument (%d)\n",
             var);
      return 0;
   }

   /* move to position in file */
   pos = grid_position( v, time, var );
   if (lseek( v->FileDesc, pos, SEEK_SET )<0) {
      /* lseek failed, return error */
      printf("Error in v5dWrite[Compressed]Grid: seek failed, disk full?\n");
      return 0;
   }

   /* write ga, gb arrays */
   k = 0;
   if (write_float4_array( v->FileDesc, ga, v->Nl[var] ) == v->Nl[var] &&
       write_float4_array( v->FileDesc, gb, v->Nl[var] ) == v->Nl[var]) {
      /* write compressed grid data (k=1=OK, k=0=Error) */
      n = v->Nr * v->Nc * v->Nl[var];
      if (v->CompressMode==1) {
         k = write_block( v->FileDesc, compdata, n, 1 )==n;
      }
      else if (v->CompressMode==2) {
         k = write_block( v->FileDesc, compdata, n, 2 )==n;
      }
      else if (v->CompressMode==4) {
         k = write_block( v->FileDesc, compdata, n, 4 )==n;
      }
   }

   if (k==0) {
      /* Error while writing */
      printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
   }
   return k;

/*
   n = v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
   if (write_bytes( v->FileDesc, compdata, n )!=n) {
      printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
      return 0;
   }
   else {
      return 1;
   }
*/
}




/*
 * Compress a grid and write it to a v5d file.
 * Input:  v - pointer to v5dstruct describing the file
 *         time, var - which timestep and variable (starting at 0)
 *         data - address of uncompressed grid data
 * Return:  1 = ok, 0 = error.
 */
int v5dWriteGrid( v5dstruct *v, int time, int var, const float data[] )
{
   float ga[MAXLEVELS], gb[MAXLEVELS];
   void *compdata;
   int n, bytes;
   float min, max;

   if (v->Mode!='w') {
      printf("Error in v5dWriteGrid: file opened for reading,");
      printf(" not writing.\n");
      return 0;
   }
   if (time<0 || time>=v->NumTimes) {
      printf("Error in v5dWriteGrid: bad timestep argument (%d)\n", time);
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Error in v5dWriteGrid: bad variable argument (%d)\n", var);
      return 0;
   }

   /* allocate compdata buffer */
   if (v->CompressMode==1) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned char);
   }
   else if (v->CompressMode==2) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned short);
   }
   else if (v->CompressMode==4) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(float);
   }
   compdata = (void *) malloc( bytes );
   if (!compdata) {
      printf("Error in v5dWriteGrid: out of memory (needed %d bytes)\n",
             bytes );
      return 0;
   }

   /* compress the grid data */
   v5dCompressGrid( v->Nr, v->Nc, v->Nl[var], v->CompressMode, data,
                    compdata, ga, gb, &min, &max );

   /* update min and max value */
   if (min<v->MinVal[var])
      v->MinVal[var] = min;
   if (max>v->MaxVal[var])
      v->MaxVal[var] = max;

   /* write the compressed grid */
   n = v5dWriteCompressedGrid( v, time, var, ga, gb, compdata );

   /* free compdata */
   free( compdata );

   return n;
}



/*
 * Close a v5d file which was opened with open_v5d_file() or
 * create_v5d_file().
 * Input: f - file descriptor
 * Return:  1 = ok, 0 = error
 */
int v5dCloseFile( v5dstruct *v )
{
   int status = 1;

   if (v->Mode=='w') {
      /* rewrite header because writing grids updates the minval and */
      /* maxval fields */
      lseek( v->FileDesc, 0, SEEK_SET );
      status = write_v5d_header( v );
      lseek( v->FileDesc, 0, SEEK_END );
      close( v->FileDesc );
   }
   else if (v->Mode=='r') {
      /* just close the file */
      close(v->FileDesc);
   }
   else {
      printf("Error in v5dCloseFile: bad v5dstruct argument\n");
      return 0;
   }
   v->FileDesc = -1;
   v->Mode = 0;
   return status;
}




/**********************************************************************/
/*****           Simple v5d file writing functions.               *****/
/**********************************************************************/



static v5dstruct *Simple = NULL;



/*
 * Create a new v5d file specifying both a map projection and vertical
 * coordinate system.  See README file for argument details.
 * Return:  1 = ok, 0 = error.
 */
int v5dCreate( const char *name, int numtimes, int numvars,
               int nr, int nc, const int nl[],
               const char varname[MAXVARS][10],
               const int timestamp[], const int datestamp[],
               int compressmode,
               int projection,
               const float proj_args[],
               int vertical,
               const float vert_args[] )
{
   int var, time, maxnl, i;

   /* initialize the v5dstruct */
   Simple = v5dNewStruct();

   Simple->NumTimes = numtimes;
   Simple->NumVars = numvars;
   Simple->Nr = nr;
   Simple->Nc = nc;
   maxnl = nl[0];
   for (var=0;var<numvars;var++) {
      if (nl[var]>maxnl) {
         maxnl = nl[var];
      }
      Simple->Nl[var] = nl[var];
      Simple->LowLev[var] = 0;
      strncpy( Simple->VarName[var], varname[var], 10 );
      Simple->VarName[var][9] = 0;
   }

   /* time and date for each timestep */
   for (time=0;time<numtimes;time++) {
      Simple->TimeStamp[time] = timestamp[time];
      Simple->DateStamp[time] = datestamp[time];
   }

   Simple->CompressMode = compressmode;

   /* Map projection and vertical coordinate system */
   Simple->Projection = projection;
   memcpy( Simple->ProjArgs, proj_args, MAXPROJARGS*sizeof(float) );

   Simple->VerticalSystem = vertical;
   if (vertical == 3) {
     /* convert pressures to heights */
     for (i=0; i<MAXVERTARGS; i++) {
       if (vert_args[i] > 0.000001) {
         Simple->VertArgs[i] = pressure_to_height(vert_args[i]);
       }
       else Simple->VertArgs[i] = 0.0;
     }
   }
   else {
     memcpy( Simple->VertArgs, vert_args, MAXVERTARGS*sizeof(float) );
   }

   /* create the file */
   if (v5dCreateFile( name, Simple )==0) {
     printf("Error in v5dCreateSimpleFile: unable to create %s\n", name );
     return 0;
   }
   else {
      return 1;
   }
}



/*
 * Create a new v5d file using minimal information.
 * Return:  1 = ok, 0 = error.  See README file for argument details.
 */
int v5dCreateSimple( const char *name, int numtimes, int numvars,
                     int nr, int nc, int nl,
                     const char varname[MAXVARS][10],
                     const int timestamp[], const int datestamp[],
                     float northlat, float latinc,
                     float westlon, float loninc,
                     float bottomhgt, float hgtinc )
{
   int nlvar[MAXVARS];
   int compressmode, projection, vertical;
   float proj_args[100], vert_args[MAXLEVELS];
   int i;

   for (i=0;i<numvars;i++) {
      nlvar[i] = nl;
   }

   compressmode = 1;

   projection = 1;
   proj_args[0] = northlat;
   proj_args[1] = westlon;
   proj_args[2] = latinc;
   proj_args[3] = loninc;

   vertical = 1;
   vert_args[0] = bottomhgt;
   vert_args[1] = hgtinc;

   return v5dCreate( name, numtimes, numvars, nr, nc, nlvar,
                     varname, timestamp, datestamp, compressmode,
                     projection, proj_args, vertical, vert_args );
}



/*
 * Set lowest levels for each variable (other than default of 0).
 * Input: lowlev - array [NumVars] of ints
 * Return:  1 = ok, 0 = error
 */
int v5dSetLowLev( int lowlev[] )
{
  int var;

  if (Simple) {
     for (var=0;var<Simple->NumVars;var++) {
        Simple->LowLev[var] = lowlev[var];
     }
     return 1;
  }
  else {
     printf("Error: must call v5dCreate before v5dSetLowLev\n");
     return 0;
  }
}


/*
 * Set the units for a variable.
 * Input:  var - a variable in [1,NumVars]
 *         units - a string
 * Return:  1 = ok, 0 = error
 */
int v5dSetUnits( int var, const char *units )
{
  if (Simple) {
     if (var>=1 && var<=Simple->NumVars) {
        strncpy( Simple->Units[var-1], units, 19 );
        Simple->Units[var-1][19] = 0;
        return 1;
     }
     else {
        printf("Error: bad variable number in v5dSetUnits\n");
        return 0;
     }
  }
  else {
     printf("Error: must call v5dCreate before v5dSetUnits\n");
     return 0;
  }
}



/*
 * Write a grid to a v5d file.
 * Input:  time - timestep in [1,NumTimes]
 *         var - timestep in [1,NumVars]
 *         data - array [nr*nc*nl] of floats
 * Return:  1 = ok, 0 = error
 */
int v5dWrite( int time, int var, const float data[] )
{
   if (Simple) {
      if (time<1 || time>Simple->NumTimes) {
         printf("Error in v5dWrite: bad timestep number: %d\n", time );
         return 0;
      }
      if (var<1 || var>Simple->NumVars) {
         printf("Error in v5dWrite: bad variable number: %d\n", var );
      }
      return v5dWriteGrid( Simple, time-1, var-1, data );
   }
   else {
      printf("Error: must call v5dCreate before v5dWrite\n");
      return 0;
   }
}



/*
 * Close a v5d file after the last grid has been written to it.
 * Return:  1 = ok, 0 = error
 */
int v5dClose( void )
{
   if (Simple) {
     int ok = v5dCloseFile( Simple );
     v5dFreeStruct( Simple );
     return ok;
   }
   else {
     printf("Error: v5dClose: no file to close\n");
     return 0;
   }
}



/**********************************************************************/
/*****                FORTRAN-callable simple output              *****/
/**********************************************************************/


/*
 * Create a v5d file.  See README file for argument descriptions.
 * Return:  1 = ok, 0 = error.
 */
#ifdef FORTRANUNDERSCORE 
   int v5dcreate_
#else
#  ifdef FORTRANDOUBLEUNDERSCORE
     int v5dcreate_
#  else
#    ifdef _CRAY
       int V5DCREATE
#    else
       int v5dcreate
#    endif
#  endif
#endif
           ( const char *name, const int *numtimes, const int *numvars,
             const int *nr, const int *nc, const int nl[],
             const char varname[][10],
             const int timestamp[], const int datestamp[],
             const int *compressmode,
             const int *projection,
             const float proj_args[],
             const int *vertical,
             const float vert_args[] )
{
   char filename[100];
   char names[MAXVARS][10];
   int i, maxnl, args;

   /* copy name to filename and remove trailing spaces if any */
   copy_string( filename, name, 100 );

   /*
    * Check for uninitialized arguments
    */
   if (*numtimes<1) {
      printf("Error: numtimes invalid\n");
      return 0;
   }
   if (*numvars<1) {
      printf("Error: numvars invalid\n");
      return 0;
   }
   if (*nr<2) {
      printf("Error: nr invalid\n");
      return 0;
   }
   if (*nc<2) {
      printf("Error: nc invalid\n");
      return 0;
   }
   maxnl = 0;
   for (i=0;i<*numvars;i++) {
      if (nl[i]<1) {
         printf("Error: nl(%d) invalid\n", i+1);
         return 0;
      }
      if (nl[i]>maxnl) {
         maxnl = nl[i];
      }
   }

   for (i=0;i<*numvars;i++) {
      if (copy_string2( names[i], varname[i], 10)==0) {
         printf("Error: unitialized varname(%d)\n", i+1);
         return 0;
      }
   }

   for (i=0;i<*numtimes;i++) {
      if (timestamp[i]<0) {
         printf("Error: times(%d) invalid\n", i+1);
         return 0;
      }
      if (datestamp[i]<0) {
         printf("Error: dates(%d) invalid\n", i+1);
         return 0;
      }
   }

   if (*compressmode != 1 && *compressmode != 2 && *compressmode != 4) {
      printf("Error: compressmode invalid\n");
      return 0;
   }

   switch (*projection) {
      case 0:
         args = 4;
         break;
      case 1:
         args = 0;
         if (IS_MISSING(proj_args[0])) {
            printf("Error: northlat (proj_args(1)) invalid\n");
            return 0;
         }
         if (IS_MISSING(proj_args[1])) {
            printf("Error: westlon (proj_args(2)) invalid\n");
            return 0;
         }
         if (IS_MISSING(proj_args[2])) {
            printf("Error: latinc (proj_args(3)) invalid\n");
            return 0;
         }
         if (IS_MISSING(proj_args[3])) {
            printf("Error: loninc (proj_args(4)) invalid\n");
            return 0;
         }
         break;
      case 2:
         args = 6;
         break;
      case 3:
         args = 5;
         break;
      case 4:
         args = 7;
         break;
/* MJK 12.12.99 */
      case 5:
         args = 4;
         break;
      default:
         args = 0;
         printf("Error: projection invalid\n");
         return 0;
   }
   for (i=0;i<args;i++) {
      if (IS_MISSING(proj_args[i])) {
         printf("Error: proj_args(%d) invalid\n", i+1);
         return 0;
      }
   }

   switch (*vertical) {
      case 0:
/* WLH 31 Oct 96  -  just fall through 
         args = 4;
         break;
*/
      case 1:
         args = 0;
         if (IS_MISSING(vert_args[0])) {
            printf("Error: bottomhgt (vert_args(1)) invalid\n");
            return 0;
         }
         if (IS_MISSING(vert_args[1])) {
            printf("Error: hgtinc (vert_args(2)) invalid\n");
            return 0;
         }
         break;
      case 2:
      case 3:
         args = maxnl;
         break;
      default:
         args = 0;
         printf("Error: vertical invalid\n");
         return 0;
   }
   for (i=0;i<args;i++) {
      if (IS_MISSING(vert_args[i])) {
         printf("Error: vert_args(%d) invalid\n", i+1);
         return 0;
      }
   }

   return v5dCreate( filename, *numtimes, *numvars, *nr, *nc, nl,
                     (const char(*)[10]) names, timestamp, datestamp,
                     *compressmode,
                     *projection, proj_args, *vertical, vert_args );
}




/*
 * Create a simple v5d file.  See README file for argument descriptions.
 * Return:  1 = ok, 0 = error.
 */
#ifdef FORTRANUNDERSCORE
  int v5dcreatesimple_
#else
#  ifdef _CRAY
     int V5DCREATESIMPLE
#  else
     int v5dcreatesimple
#  endif
#endif
           ( const char *name, const int *numtimes, const int *numvars,
             const int *nr, const int *nc, const int *nl,
             const char varname[][10],
             const int timestamp[], const int datestamp[],
             const float *northlat, const float *latinc,
             const float *westlon, const float *loninc,
             const float *bottomhgt, const float *hgtinc )
{
   int compressmode, projection, vertical;
   float projarg[100], vertarg[MAXLEVELS];
   int varnl[MAXVARS];
   int i;

   for (i=0;i<MAXVARS;i++) {
      varnl[i] = *nl;
   }

   compressmode = 1;

   projection = 1;
   projarg[0] = *northlat;
   projarg[1] = *westlon;
   projarg[2] = *latinc;
   projarg[3] = *loninc;

   vertical = 1;
   vertarg[0] = *bottomhgt;
   vertarg[1] = *hgtinc;

#ifdef FORTRANUNDERSCORE
   return v5dcreate_
#else
#  ifdef FORTRANDOUBLEUNDERSCORE
     return v5dcreate_
#  else
#    ifdef _CRAY
       return V5DCREATE
#    else

       return v5dcreate
#    endif
#  endif
#endif
                   ( name, numtimes, numvars, nr, nc, varnl,
                     varname, timestamp, datestamp, &compressmode,
                     &projection, projarg, &vertical, vertarg );
}



/*
 * Set lowest levels for each variable (other than default of 0).
 * Input: lowlev - array [NumVars] of ints
 * Return:  1 = ok, 0 = error
 */
#ifdef FORTRANUNDERSCORE
   int v5dsetlowlev_
#else
#  ifdef _CRAY
     int V5DSETLOWLEV
#  else
     int v5dsetlowlev
#  endif
#endif
          ( int *lowlev )
{
   return v5dSetLowLev(lowlev);
}



/*
 * Set the units for a variable.
 * Input: var - variable number in [1,NumVars]
 *        units - a character string
 * Return:  1 = ok, 0 = error
 */
#ifdef FORTRANUNDERSCORE
   int v5dsetunits_
#else
#  ifdef _CRAY
     int V5DSETUNITS
#  else
     int v5dsetunits
#  endif
#endif
          ( int *var, char *name )
{
   return v5dSetUnits( *var, name );
}



/*
 * Write a grid of data to the file.
 * Input:  time - timestep in [1,NumTimes]
 *         var - timestep in [1,NumVars]
 *         data - array [nr*nc*nl] of floats
 * Return:  1 = ok, 0 = error
 */
#ifdef FORTRANUNDERSCORE
   int v5dwrite_
#else
#  ifdef FORTRANDOUBLEUNDERSCORE
     int v5dwrite_
#  else
#    ifdef _CRAY
       int V5DWRITE
#    else
       int v5dwrite
#    endif
#  endif
#endif
          ( const int *time, const int *var, const float *data )
{
   return v5dWrite( *time, *var, data );
}



/*
 * Specify the McIDAS GR3D file number and grid number which correspond
 * to the grid specified by time and var.
 * Input:  time, var - timestep and variable of grid (starting at 1)
 *         mcfile, mcgrid - McIDAS grid file number and grid number
 * Return:  1 = ok, 0 = errror (bad time or var)
 */
#ifdef FORTRANUNDERSCORE
   int v5dmcfile_
#else
#  ifdef _CRAY
     int V5DMCFILE
#  else
     int v5dmcfile
#  endif
#endif
         ( const int *time, const int *var,
           const int *mcfile, const int *mcgrid )
{
   if (*time<1 || *time>Simple->NumTimes) {
      printf("Bad time argument to v5dSetMcIDASgrid: %d\n", *time );
      return 0;
   }
   if (*var<1 || *var>Simple->NumVars) {
      printf("Bad var argument to v5dSetMcIDASgrid: %d\n", *var );
      return 0;
   }

   Simple->McFile[*time-1][*var-1] = (short) *mcfile;
   Simple->McGrid[*time-1][*var-1] = (short) *mcgrid;
   return 1;
}



/*
 * Close a simple v5d file.
 */
#ifdef FORTRANUNDERSCORE
   int v5dclose_( void )
#else
#  ifdef FORTRANDOUBLEUNDERSCORE
     int v5dclose_( void )
#  else
#    ifdef _CRAY
       int V5DCLOSE( void )
#    else
       int v5dclose( void )
#    endif
#  endif
#endif
{
   return v5dClose();
}
